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Abstract 

We use dense evaluation ordering to define HRD (Hybrid- Restriction Diagram), a new BDD-likc 
data-structure for the representation and manipulation of state-spaces of linear hybrid automata. Wc 
present and discuss various manipulation algorithms for HRD. including the basic set-oriented opera- 
tions, weakest precondition calculation, and normalization. Wc implemented the ideas and experimented 
to see their performance. Finally, we have also developed a pruning technique for state-space explo- 
ration based on parameter valuation space characterization. The technique showed good promise in our 
experiment. 
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1 Introduction 

Symbolic analysis of linear hybrid automata (LHA) [5, 7] can generate a symbolic characterization of the 
reachable state-space of the LHA. When static parameters (system variables whose values are decided before 
| run-time and never changed in run-time) are used in LHAs, such symbolic characterizations may shed 
important feedback information to engineers. For example, we may use such symbolic characterizations 
to choose proper parameter values to avoid from unsafe system designs. Unfortunately, LHA systems are 
extremely complex and not subject to algorithmic analysis [8]. Thus in real- world applications, it is very 
important to use every measure to enhance the efficiency of LHA parametric analysis. 

In this work, we extend BDD-like data-structures [10, 12] for the representation and manipulation of 
LHA state-spaces. BDD-like data-structures have the advantage of data-sharing in both representation 
and manipulation and have shown great success in VLSI verification industry. One of the major difficulties 
to use BDD-like data-structures to analyze LHAs comes from the unboundedness of the dense variable 
value ranges and the unboundedness of linear constraints. To explain one of the major contribution of 
this work, we need to discuss the following issue first. In the research of BDD-like data-structures, there 
are two classes of variables: system variables and decision atoms [23]. System variables are those used 
in the input behavior descriptions. Decision atoms are those labeled on each BDD nodes. For discrete 
systems, these two classes are the same, that is, decision atoms are exactly those system variables. But 
for dense-time systems, decision atoms can be different from state varaibles. For example, in CDD [11] 
and CRD [23], decision atoms are of the form x — x' where x and x' are system variables of type clock. 
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Previous work on BDD-like data-structures are based on the assumption that decision atom domains are of 
finite sizes. Thus we need new techniques to extend BDD-like data-structures to represent and manipulate 
state-spaces of LHAs. Our innovations include using constraints like — 3A + x — Ay (where A, x, y are 
dense variables), as the decision atoms and using total dense orderings among these atoms. In this way, 
we devised HRD (Hybrid-Restriction Diagram) and successfully extend BDD-technology to models with 
unbounded domains of decision atoms. 

In total, we defined three total dense-orderings for HRD constriants (section 6). We also present 
algorithms for set-oriented operations (section 7) and symbolic weakest precondition calculation (section 9), 
procedures for symbolic parametric analysis (section 9), and discuss our implementation of symbolic convex 
polyhedra representation normalization (section 10). Especially, in the presentation of our previous work 
of BDD-like data-structures for timed automata, people usually asked for presentation of our algorithms 
for weakest precondition construction. In this paper, we endeavored to make a concise presentation. 

We have also developed a techique for fast parametric analysis of LHA (section 11). The technique 
prunes state-space exploration based on static parameter space characterization. The technique gives us 
very good performance. Desirably, this technique does not sacrifice the precision of parametric analysis. 
Especially, for one benchmark, the state-space exploration does not converge without this technique! To 
our knowledge, nobody else has come up with a similar technique. Finally, we have implemented our ideas 
in our tool red 5.0 and reported our experiments to see how the three dense-orderings perform and how 
our implementation performs in comparison with HyTech 2.4.5 [14] and TReX 1.3 [1,3]. 

2 Related work 

Many modern model-checkers [18, 23, 27] for timed automata [6] are built around symbolic manipulation 
procedures [7, 15] of zones, which means behaviorally equivalent convex state spaces of timed automata. 
The most popular data-structure for zones is DBM [13], which is a two dimensional matrix recording 
differences between pairs of clocks and nothing BDD-like. 

As far as we know, the first paper that discusses how to use BDD to encode zones is by Wang, Mok, 
and Emerson in 1993 [25]. They discussed how to use BDD with decision atoms like Xi + c < Xj + d to 
model-check timed automata. Here c and d are timing constants with magnitude < Ca- However, they 
did not report implementation and experiments. In the last several years, people have explored in this 
approach in the hope to duplicate the success of BDD techniques [10, 12] in hardware verification for the 
verification of timed automata [2, 9, 11, 16, 17, 19-23]. 

For parametric analysis, Annichini et al have extended DBM to PDBM for parametric analysis of timed 
automata [1,3] and implemented a tool called TReX, which also supports verification with lossy channels. 
Due to the differences in their target systems, it can be difficult to directly compare the performances of 
TReX and our implementation red 5.0. For example, TReX only allows for clocks while red 5.0 allows 
for dense variables with rate intervals. To construct time-progress weakest preconditions (or strongest 
postcondition in forward analysis) for systems with dense variable rate intervals, red 5.0 needs to use 
one (5-variable for each dense variables and significantly increase the number of decision atoms involving 
5- variables. According to the new formulation of time-prorgress weakest precondition algorithm in [23], for 
systems with only clocks, no literals involving <5-variables ever need to be generated. Thus the complexity 
for the algorithm used in red 5.0 is relatively higher than those used in TReX. On the other hand, TReX 
may have tuned its performance for the verification of lossy channel systems. 

For LHAs, people also used convex subspaces, called convex polyhedra, as basic unit for symbolic manip- 
ulation. A convex polyhedron characterizes a state-space of an LHA and can be symbolically represented 
by a set of constraints like a\X\ + . . . + a n x n ~ c [4, 5, 7]. Two commonly used representations for convex 
polyhedra in HyTech are (1) polyhedras and (2) frames in dense state-space [14]. These two representa- 
tions neither are BDD-like nor can represent concave state-spaces. Data-sharing among convex polyhedra 
is difficult. 
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Figure 1: Fischer's timed mutual exclusion algorithm in LHA 

3 Parametric analysis of linear hybrid automata (LHA) 

A linear hybrid automata (LHA ) [7] is a finite-state automaton equipped with a finite set of dense variables 
which can hold real- values. At any moment, the LHA can stay in only one mode (or control location). In its 
operation, one of the transitions can be triggered when the corresponding triggering condition is satisfied. 
Upon being triggered, the LHA instantaneously transits from one mode to another and sets some dense 
variables to values in certain ranges. In between transitions, all dense variables increase their readings at 
rates determined by the current mode. 

For convenience, given a set Q of modes and a set X of dense variables, we use P(Q, X) as the set of 
all Boolean combinations of atoms of the forms q and ^ a^Xi ~ c, where q G Q, a* are integers constants, 
ii G X, "~" is one of <,<,=,>, >, and c is a rational constant. 

We also let be the set of rational intervals like (d, d') where '(' is either '[' or '('; '}' is either ']' or ')'; 
and d,d' are — oo, oo, or rational numbers. 

Definition 1 linear hybrid automata (LHA) An LHA A is a tuple 
(X, Q, I, /i, 7, E, t, tt) with the following restrictions. X is the set of dense variables. Q is the set of 
modes. I G P(Q, X) is the initial condition, fj, : Q i— ► P(0, X) defines the invariance condition of each 
mode. 7 : (Q x X) \— > defines the rate intervals of dense variables. E C Q x Q is the set of transitions. 
t : E i— > P(0, X) defines the triggering condition of transitions, it is a partial function from E x X to 
that defines the interval assignments to dense variables at each transition. If ir(e,x) is undefined, x is 
not assigned a value in transition e; otherwise, x is nondeterministically assigned a finite rational value in 
7r(e, x) in transition e. ■ 

In figure 1, we have drawn a version of the Fishcer's mutual exclusion algorithm for a process. There 
are two static parameters a and (3 that controls the behavior of the processes. In each mode, local clock 
x increases its reading according to a rate in [4/5, 1]. The rate interval in each mode can be different. 

A valuation of a set is a mapping from the set to another set. Given an r] G P(Q, X) and a valuation v 
of X, we say v satisfies r/, in symbols v (= r], iff it is the case that when the variables in r\ are interpreted 
according to u, rj will be evaluated true. 

Definition 2 states A state v of A = {X, Q, I, /x, 7, 12, r, 7r) is a valuation of A U Q s.t. 

• there is a unique q £ Q such that z/(g) = true and for all q' 7^ v(q') = /aZse; 

• for each x £ X, v{x) 6 7£ (the set of reals) and \/q € Q, f (<?) =>■ v \= j-t(q). 

Given state v and q 6 Q such that i/(g) = irue, we call g the mode of z^, in symbols 1/^. ■ 
For any t € 1Z + (the set of nonnegative reals), v -w 1/ iff we can go from to z/ merely by the passage 

of t time units. Formally speaking, v -w z/ is true iff 1/ is a state identical to z^ except that for every i£l 
with 7(z/ Q , 2) = (d, d'), v'(x) G (z/(x) + 1 ■ d, v(x) + t-d'). 

For a transition e G E, v z/ iff we can go from z/ to v 1 with discrete transition e = ((7, </). Formally 
speaking, v v' is true iff 1^ = g, v \= fi(q) A r(e), and v' is identical to v except that 

• p'Q = q' and z/ [= fi(q'); and 

• for each x E X, if 7r(e,rE) is defined, z/(x) g vr(e, sc); otherwise, z/'(x) = z^(x); 

Definition 3 runs Given an LHA A = {X, Q, I, fi, 7, E, r, tt), a run is an infinite sequence of pairs 

(</>o, to)((j)i, ti) . . . (4>k,tk) such that £0*1 • • -tk is a monotonically increasing real-number (time) 

divergent sequence and for all k > 0, 

• cj)k is a mapping from [tk,tk+i] to states, and 
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• time-progress is continuous: that is, 

t'—t 

for each t\. < t < t' < t k +i, 4>k{t) ~* 4>k{t')', and 

• invariance conditions are preserved in each interval: that is, 
for all t k <t < tk+i, 4>k(t) h Li(4>k(t) Q ); and 

• either no transition happens at time tfc+i, that is, (pk(tk+i)® = (pk+iitk+i)®', or a transition e happens 
at t k+ i, that is, <j) k (tk+i) ^ 4>k+i{tk+i) ■ ■ 

A run p = (cfto, to)(4>i,ti) . . . {4>k,tk) is safe w.r.t. a safety state-predicate 77, in symbols p \= 77, iff 

for all k > and t € [t k , tfc+i], <j>k{t) \= A dense variable x in an LHA is a static parameter iff its rate is 
always zero in all modes. Suppose H is the set of static parameters in X of LHA A. A static parameter 

valuation TC of a run (</>o, toX^i) *i) • • • {4>k-> tk) is a mapping from H to reals such that TC is consistent 

with every state along p, i.e., Vx € -ffV/c > 0(4>k(tk)( x ) = TL{x)). TC is a parametric solution to ^4 and 77 iff 
for all runs /) with static parameter valuation 7i, p \= n. 

Our verification framework is called parametric safety analysis problem. A parametric safety analysis 
problem instance, PSA(j4, 77) in notations, consists of an LHA A and a safety state-predicate 77 E P(Q, X). 
Such a problem instance asks for a symbolic characterization of all parametric solutions to A and 77. The 
general parametric safety analysis problem is undecidable. 

4 Convex polyhedra 

Given a set X = {xi, . . . ,x n } of dense system variables, an LH-expression (linear hybrid expression) is 
an expression like a\X\ + . . . + a n x n where a\,...,a n are integer constants. It is normalized iff the gcd 
of nonzero coefficients in {a±, . . . ,a n } is 1, i.e., gcd{oj | 1 < i < n; a, 7^ 0} = 1. From now on, we shall 
assume that all given LH-expressions are normalized. 

An LH-upperbound is either (<,oo) or a pair like (~,c) where ~€ {" < "," < "} and c is a rational 
number. There is a natural ordering C among the LH-upperbounds. That is for any two (~, c) and (~', c'), 
(~,c) C (~',c') iff c < d or (c = c'A ~= " < "A ~'= " < "). Intuitively, if (~,c) C (~',c'), then (~,c) is 
more restrictive than (~',c'). 

An LH-constraint is a pair of an LH-expression and an LH-upperbound. Given an LH-expression 
ajXj and an LH-upperbound (~ c), we shall naturally write the corresponding LH-constraint as 
aiXi ~ c. A convex polyhedron is symbolically represented by a conjunction of LH-constraints and 
means a behaviorally equivalent state subspace of an LHA. Formally, a convex polyhedron £ can be de- 
fined as a mapping from the set of LH-expressions to the set of LH-upperbounds. Alternatively, we may 
also represent a convex polyhedron ( as the set aiXi ~ c | CEi a i x i) = (~j c )}- We shall use the two 
equivalent notations flexibly as we see fit. With respect to a given X, the set of all LH-expressions and 
the set of convex polyhedra are both infinite. 

5 HRD (Hybrid-Restriction Diagram) 

To construct BDD-like data-structures, three fundamental issues have to be solved. The first is the 
domain of the decision atoms; the second is the range of the arc labels from BDD nodes; and the third 
is the evaluation ordering among the decision atoms. For modularity of presentation, we shall leave the 
discussion of the evaluation orderings to section 6. In this section, we shall assume that we are given a 
decision atom evaluation ordering. 

We decide to follow an approach similar to the one adopted in [23]. That is, our decision atoms 
will be LH-expressions while our BDD arcs will be labeled with LH-upperbounds. A node labeled with 
decision atom £V aixi together with a corresponding outgoing arc label (~, c) constitute the LH-constraint 
of Y2i a i x i ~ c. A root-to-terminal path in an HRD thus represents the conjunction of constituent LH- 
constraints along the path. Figure 2(a) is an example of our proposed BDD-like data-strucutre for the 
concave space of 

(x 2 -x 3 < -5/7 V -5A - 2x 2 + 10x 3 < 48/7) A A - x 2 + 10x 3 < 9 
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Figure 2: Examples of HRD 



assuming that —5^4 — 2x2 + IOX3 precedes X2 — X3 (in symbols —5A — 2x2 + IOX3 -< X2 — X3) and X2 — X3 
precedes A — X2 + IOX3 in the given evaluation ordering. In this example, the system variables are A, X2, £3 
while the decision atoms are X2 — £3, —5^4 — 2x2 + IOX3, and A — X2 + IOX3. 

Definition 4 HRD (Hybrid- Restriction Diagram) Given dense variable set X = {xi, . . . , x n } and an evalu- 
ation ordering ~< among normalized LH-expressions of X , an HRD is either true or a tuple (v, (f3\ , D\ ) , . . . , (/3 m , D m ) ) 
such that 

• v is a normalized LH-expression; 

• for each 1 < i < m, Pi is an LH-upperbound s.t. (<, 00) 7^ j3\ C P2 C • • • C j9 m ; and 

• for each 1 < i < m, -Dj is an HRD such that if Dj = (i>j, . . .), then ?; -< Uj. 

For completeness, we use "()" to represent the HRD for false. ■ 
In our algorithms, false does not participate in comparison of evaluation orderings among decision 
atoms. Also, note that in figure 2, for each arc label (~,c), we simply put down ~ c for convenience. 
Note that an HRD records a set of convex polyhedra and each root-leaf path represents such a convex 
polyhedron. 

6 Three dense orderings among decision atoms 

In the definition of a dense-ordering among decision atoms (i.e., LH-expressions), special care must be taken 
to facilitate efficient manipulation of HRDs. Here we use the experience reported in [23] and present three 
criteria in designing the orderings among LH-expressions. The three criteria are presented in sequence 
proportional to their respective importances. 

First, it is desirable to place a pair of converse LH-expressions next to one another so that simple 
inconsistencies can be easily detected. That is, LH-expressions Yli a i x i an d Yli ~ a i x i are better placed 
next to one another in the ordering. For example, with this arrangement, the inconsistency of —x\ +3x2 < 
—5 Aii - 3x2 < can be checked by comparing adjacent nodes in HRD paths. To fulfill this requirement, 
when comparing the precedence between LH-expressions in a given ordering, we shall first toggle the signs 
of coefficients of an LH-expression if its first nonzero coefficient is positive. If two LH-expressions are 
identical after the necessary toggling, then we compare the signs of their first nonzero coefficients to decide 
the precedence between the two. 

With the requirement mentioned in the last paragraph, from now on, we shall only focus on the 
orderings among LH-expressions whose first nonzero coefficients are negative. 

Secondly, according to past experience reported in the literature, it is important to place strongly 
correlated LH-expressions close together in the evaluation orderings. Usually, instead of a single global 
LHA, we are given a set of communicating LHAs, each representing a process. Thus it is desirable to 
place LH-expressions for the same process close to each other in the orderings. Our second important 
criterion respects this experience. Given a system with m processes with respective local dense variables, 
we shall partition the LH-expressions into m + 1 groups: Gq,G\, . . . , G m . Gq contains all LH-expressions 
without local variables (i.e., coefficients for local variables are all zero). For each p > 0, G p contains 
all LH-expressions with a nonzero coefficient for a local variable of process p and only zero coefficients 
for local variables of processes p + 1, . . . ,m. Then our second criterion requires that for all < p < m, 
LH-expressions in G p precede those in G p+ ±, . . . , G m . 

If the precedence between two LH-expressions cannot be determined with the two above-mentioned cri- 
teria, then the following third criterion comes to play. This one is a challenge since for each of Go, ... , G m 
can be of infinite size. Traditionally, BDD-like data-structures have been used with finite decision atom 
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domains. But now we need to find a way to determine the precedence among infinite number of LH- 
expressions (our decision atoms in HRD). For this purpose, we invent to use dense-orderings among LH- 
expressions. We shall present three such orderings in the following. Sometimes it is difficult to predict 
which orderings are better suitable for what kind of verification tasks. In section 12, we shall report 
experiments with these orderings. 

Dictionary ordering: We can represent each LH-expression as a string, assuming that the ordering 
among x\, . . . , x n is fixed and no blanks are used in the string. Then we can use dictionary ordering and 
ASCII ordering to decide the precedence among LH-expressions. For the LH-expressions in figure 2, we 
then have — 5A — 2x2 + IOX3 -< A — xi + 10x3 -< X2 — £3 since '— ' precedes 7 A' and 'A' precedes 'x' in 
ASCII. The corresponding HRD in dictionary ordering is in figure 2(c). One interesting feature of this 
ordering is that it has the potential to be extended to nonlinear hybrid constraints. For example, we may 
say cos(xi) + x| ■< x\ — X2X3 in dictionary ordering since 'c' precedes 'x' in ASCII. 

Coefficient ordering: Assume that the ordering of the dense variables is fixed as x\, . . . ,x m . In this 
ordering, the precedence between two LH-expressions is determined by iteratively comparing the co- 
efficients of dense variables x\,...,x n in sequence. For the LH-expressions in figure 2, we then have 
—5A — 2x2 + IOX3 ~< X2 — X3 ~< A — X2 + IOX3 The HRD in this ordering is in figure 2(a). 

Magnitude ordering: This ordering is similar to the last one. Instead of comparing coefficients, we 
compare the absolute values of coefficients. We iteratively 

• first compare the absolute values of coefficients of Xj, and 

• if they are equal, then compare the signs of coefficients of Xj. 

For the LH-expressions in figure 2, we then have X2 — X3 -< A — X2 + IOX3 ~< —5A — 2x2 + IOX3 in this 
magnitude ordering. The HRD in this ordering is in figure 2(b). 

7 Set-oriented operations 

Please be reminded that an HRD records a set of convex polyhedra. For convenience of discussion, given an 
HRD, we may just represent it as the set of covnex polyhedra recorded in it. Definitions of set-union (U), 
set-intersection (n), and set-exclusion (— ) of two convex polyhedra sets respectively represented by two 
HRDs are straightforward. For example, given HRDs D\ : {(j, £2} and D2 : {C2, C3 }: D\ n D2 is the HRD 
for {C2}; Di U D2 is for {Ci>C2jC3}; an d D\ — D2 is for {Ci}- The complexities of the three manipulations 
are all 0(\Di\ ■ \D 2 \). 

Given two convex polyhedra £1 and £2; (jFl £2 is a new convex polyhdron representing the space- 
intersection of Ci and £2- Formally speaking, for decision atom Yli a i x ii Ci n C2(^i a « x «) = Ci(Sj a « x «) ^ 
Ci(X^ a « x «) E C2(X^ a « x «); or C2^2,iCLiXi) otherwise. Space-intersection (n) of two HRDs D\ and D2, in 
symbols D\ n D 2 , is a new HRD for {Ci n C2 I Ci € A.; C2 € D 2 }. 

Given an evaluation ordering, we can write HRD-manipulation algorithms pretty much as usual [10,12, 
19,23]. For convenience of presentation, we may repersent an HRD (u, (@i,Bi), . . . , (f3 n , B n )) symbolically 
as (u, (/3i,Bi)i<i< n ). A union operation L)(B,D) can then be implemented as follows. 

set /* database for the recording of already-processed cases */ 
U(B,D) { 

if B = false, return D; else if D = false, return B; 
^ := 0; return recU(B,D); 

} 

rec\j(B,D) where B = (u, #i)i<i< n ), D = (v, (atj, Dj)x<j< m ) { 

if B is true or D is true, return true; else if 3F, (B, D, F) £ return F; (1) 

else if u -<, v, construct F := (u, (Pi, recU(£?j, -D))i<i<n); 
else if v -< u, construct F := (v, (aj,recU(B,Dj))i<j< m ); 
else { 

% := n;j := m; F := false; 
while i > 1 and j > 1, do { 
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if Pi = aj, { 



B := B U B { ; D : 

} 



D U Dj- F:=FU («, (ft, recU(B i( * 



— ;j — ; 



else if Pi \Z (Xj, { F 
else if ctj \Z Pi, { F 



FU(u, ( aj ,Dj));j :=j-l; } 
FU(u, (p i ,B i ));i:=i-l; } 



} 

ifi>l, F:=FU(u, (/9i,Bi)i< fc <i); 
if j > 1, F -Fuju^ai^i)^-); 

} 

* := * U {{B, D, F)}; return F; 



(2) 



} 



Note that in statement (1), we take advantage of the data-sharing capability of HRDs so that we do 
not process the same substructure twice. The set of ^ is maintained in statement (2). The algorithms 
for n and — are pretty much the same. The one for space intersection is much more involved and is not 
discussed here due to page- limit. 



As reported in the experiment with CRD (Clock- Restriction Diagram) [23], significant performance im- 
provement can be obtained if an integrated BDD-like data-structure for both dense constraints and discrete 
constraints is used instead of separate data-structure for them. It is also possible to combine HRD and 
BDD into one data-structure for fully symbolic manipulation. Since HRD only has one sink node: true, it 
is more compatible with BDD without FALSE terminal node which is more space-efficient than ordinary 
BDD. There are two things we need to take care of in this combination. The first is about the interpreta- 
tion of default values of decision atoms. In BDD, when we find a decision atom is missing during valuating 
variables along a path, the atom's value can be interpreted as either TRUE or FALSE. But in HRD, when 
we find a decision atom ^ i cnxi is missing along a path, then the constraint is interpreted as ^ cnxi < oo. 

The second is about the interpretation of HRD manipulations to BDD decision atoms. Straightfor- 
wardly, "U" and "H" on BDD decision atoms are respectively interpreted as "V" and "A" on BDD decision 
atoms. Di — D2 on BDD decision atoms is interpreted as D\ A -1-D2 when the root variable of either D\ or 
Z?2 is Boolean. For D\ n D2, the manipulation acts as "A" when either of the root are labeled with BDD 
decision atoms. Due to page-limit, we shall omit the proof for the soundness of such an interpretation. 
From now on, we shall call it HRD+BDD a combination structure of HRD and BDD. 

Finally, it is also important to define the evaluation orderings between BDD decision atoms and HRD 
decision atoms. Due to page-limit, we shall adopt the wisdom reported in [23] and place BDD decision 
atoms and HRD decision atoms that are strongly related to the same process close to each other. 

9 Weakest preconditon calculation and symbolic parametric analysis 

Our tool red runs in backward reachability analysis by default. Due to page-limit, we shall only present the 
algorithm in symbolic fashion without details. Suppose we are given an LHA A = (X, Q, I, fj,, 7, E, t, it). 
There are two basic procedures in this analysis procedure. The first, xtion(Z), e), computes the weakest 
precondition from state-space represented by HRD D through discrete transition e = (q, q'). Assume that 
the dense variables that get assigned in e are yi, ■ ■ ■ ,Vk and there is no variable that gets assigned twice 
in e. The characterization of xtion(D,e) is 



Assume that delta_exp(D) is the same as D except that all dense variables x are replaced by 
x + 5 X respectively. Here 5 X represents the value-change of variable x in time-passage. For example, 
delta_exp(2xi — 3^2 < 3/5) = 2x± + 25 X1 — 3x2 — 35 X2 < 3/5. Intuitively, when x represents the value of 

1 ye[d,d']=d<y<d'. y G {d,d'\ = d < y < d! . y G [d, d') = d < y < d' . y G (d,d') = d < y < d'. 
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H(q) n r(e) n 3yi . . . 3y k {D n ni<i< fc y,; G vr(e, yi)) 
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set R, £>; 

xtivity(Z), x) { i? := 0; return rec_xtivity(D); } 
rec_xtivity(D) { 

if D is true or false, return D; else if 3(D, D') £ R, return D'; 
else /* assume D = (ax + e, (ft, Di), ... , (/3 m ,D m )) */ { 

S : = 0; D' := Ui<i< m ax + e A n rec_xtivity_given(A, ax + e, ft); 
R:=RU {(£>, D 7 )}; return D'; 

} 

} 

rec_xtivity^given(L>, ax + e, /3) { 

if Z) is irae or /a/se, return D; else if 3(D, D') € S 1 , return D'; (3) 

else /* assume D = (bx + e', (ft, D{), (/3 m , An)) */ { 
if ab < 0, 

, _ | | / bx + e'/3j n rec_xtivity_given(L>j, ax + e, (3) \ 

:=Ul ^ m V n |6|e/gcd(o,6) + |o|e'/gcd(a,6)((|6| / 9 + |a| / a i )/gcd(a,6)) )' 
else D' := Ui<i<m ^ x + e 'A n rec_xtivity_given(Dj, ax + e, /?); 

5 := S U {(D,"D')}; return D'; (4) 

} 

(\b\(3 + [a|ft)7 § C( i(a, 5) is a shorthand for the new upperbound obtained from the xtivity of ax + e(3 and 
te + e'ft. 

Table 1: Algorithm for xtivity () 



variable x in the weakest precondition of time passage, then x + 5 X is the value of x in the postcondition 
of the time-passage. 

The second basic procedure, t±me(D,q), computes the weakest precondition from D through time 
passage in mode q. It is characterized as 

5 > n delta_exp(L>) \ 

n n i<i<n;7( g ,x 4 )=(di,^)^i € (^6,(1^5) J 

One basic building block of both xtionQ and timeQ is for the evaluation of 3x(D(x)). We implement 
this basic operation with the following symbolic procedure. 

3x(D(x)) = var_del(xtivity(L>, x), {x}). 

Procedure var_del(D, X') eliminates all constraints in D involving variables in set X' . Procedure xtivity(Z), x) 
adds to a path every constraint that can be transitively deduced from two peer constraints involving x in 
the same path in D. The algorithm of xtivityQ is in table 1. Thus we preserve all constraints transi- 
tively deducible from a dense variable before it is eliminated from a predicate. This guarantees that no 
information will be unintentionally lost after the variable elimination. 

Note that in our algorithm, we do not enumerate all paths in HRD to carry out this least fixpoint 
evaluation. Instead, in statement (3), our algorithm follows the traditional BDD programming style which 
takes advantage of the data-sharing capability of BDD-like data-structures. Thus our algorithm does 
not explode due to the combinatorial complexity of path counts in HRD. This can be justified by the 
performance of our implementation reported in section 12. 

Assume that the unsafe state is in mode qf. With the two basic procedures, then the backward 
reachable state-space from the risk state —>rf (represented as an HRD) can be characterized by 

lfpZ. (time (-177, Qf) u [Je=(q,g')eE time(xtion(Z, e), q)) 

Here If pZ.F(Z) is the least fixpoint of function FQ and is very commonly used in the reachable state- 
space representation of discrete and dense-time systems. After the fixpoint is successfully constructed, 
we conjunct it with the initial condition and then eliminate all variables except those static parameters 
(formally speaking, projecting the reachable state-space representations to the dimensions of the static 
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parameters). Suppose the set of static dense parameters is H. The characterization of unsafe parameter 
valuatons is thus 

var_del(7 n IfipZ. (time (-177, qf) U U e =(g q')eE time(xtion(Z, e), q)), X — H) 
The set of parametric solutions is characterized by the complement of this final result. 

10 Normalization 

There can be infinitely many LH-constraint sets that represent a given convex polyhedron. An LH- 
constraint in such a representation can also be redundant in that a no less restrictive upperbound can 
be derived for its LH-expression from peer LH-constraints in the same representation. To control the 
redundancy caused by recording many LH-constraint sets for the same convex polyhedron, representations 
of convex polyhedra have to be normalized. Due to page-limit, we shall skip much details in this regard. 
We emphasize that much of our implementation effort has been spent in this regard. We use a two-phase 
normalization procedure in each iteration of the least fixpoint evaluation. 

Step I, for subsumed polyhedra elimination : This step eliminates those convex polyhedra con- 
tained by a peer convex polyhedron in the HRD for the reachable state-space. First, we collect the 
LH-expressions that occur in the current reachable state-space HRD and call them proof- obligations. 
Then we try to derive the tightest constraints for these proof-obligations along each HRD paths of 
the reachable state-space representation. Then we eliminate those paths which is subsumed by other 
paths. The subsumption can be determined by pairwise comparison of all LH-constraints along two 
paths. 

Step II, for redundant constraint elimination : Along each path, we combinatorially use up to four 
constraints to check for the redundancy of peer constraints in the same path and eliminate them if 
they are found redundant. 

Again, our algorithm does not enumerate paths in HRD. Instead, it takes advantage of data-sharing 
capability of HRD for efficient processing. 

11 Pruning strategy based on parameter space construction (PSPSC) 



We have also experimented with techniques to improve the efficiency of parametric analysis. One such 
technique, called PSPSC, is avoiding new state-space exploration if the exploration does not contribute 
to new parametric solutions. A constraint is static iff all its dense variables are static parameters. Static 
constraints do not change their truth values. Once a static constraint is derived in a convex polyhedron, 
its truth value will be honored in all weakest preconditions derived from this convex polyhedron. All 
states backwardly reachable from a convex polyhedron must also satisfy the static constraints required in 
the polyhedron. Thus if we know that static parameter valuation 7i is already in the parametric solution 
space, then we really do not need to explore those states whose parameter valuations fall in H. 

With PSPSC, our new parametric analysis procedure is shown in table 2. In the procedure, we use 
varaible P to symbolically accumulate the parametric evaluations leading to the risk states in the least 
fixpoint iterations. In statement (5), we check and eliminate in D those state descriptions which cannot 
possibly contribute to new parametric evaluations by conjuncting D with —>P. 

One nice feature of PSPSC is that it does not sacrifice the precision of our parametric analysis. 

Lemma 1 TL is a parametric solution to A and n iffTi satisfies the return result o/PSA_with_PSPSC(^4, 77). 
Proof : Details omitted due to page-limit. The basic idea is that the intersection at line (5) in table 2 
only stops the further exploration of those states that do not contribute to new parameter-spaces. Those 
parameter-spaces pruned in line (5) do not contribute because they are already contained in the known 
parameter constraints P and along each exploration path, the parameter constraints only get restricter. 
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PSA_with_PSPSC(^, rj) { 

D := time(-iry, qf); D := false; P := var_del(D, X — H); 
while D ^ false, do { 
D:= DUD; 

D '■= Ue=(q,q>)eE time ( xti0Ii (.D,e),q); 

5:=5n(nP)n(nD); (5) 

P := P U var_del(I n D, X - H); 

} 

return —>P; 

1 

Table 2: Procedure for parametric safety analysis with PSPSC 



As mentioned in the proof sketch, PSPSC can help in pruning the space of exploration in big chunks. 
But in the worst case, PSPSC does not guarantee the exploration will terminate. In section 12, we shall 
report the performance of this technique. Especially, for one benchmark, the state-space exploration cannot 
converge without PSPSC. 

12 Implementation and experiments 

We have implemented our ideas in our tool red which has been previously reported in [19-23] for the 
verification of timed automata based on BDD-like data-structures, red version 5.0 supports full TCTL 
model-checking/simulation with graphical user-interface. Coverage estimation techniques for dense-time 
state-spaces has also been reported [24]. 

12.1 Comparison with HyTech 2.4.5 

We have also carried out experiments to compare various ideas mentioned in this work. In addition, 
we have also compared with HyTech 2.4.5 [14], which is the best known and most popular tool for the 
verification of LHA due to its pioneering importance. The following three benchmark series are all adapted 
from HyTech benchmark repository. 

• Fischer's mutual exclusion algorithm. This is one of the classic benchmarks. There are two static 
parameters A and B, m processes, and one local clock for each process. The first process has a local 
clock with rate in [4/5,1] while all other processes have local clocks with rates in [1,11/10]. The 
algorithm may violate the mutual exclusion property when —A < A — 11 A + 8B < 0. 

• General railroad crossing benchmarks. There is a static parameter CUTOFF, a gate-process, a controller- 
process, and m train-processes. The local dense variable of the gate-process models the angle of the 
gate and has rates in [0,0], [—10, —9], and [9, 10] depending on which modes the gate-process is in. 
The controller process does not use clocks. Each train-process uses a local clock with rate in [1, 1]. 
The system may not lower the gate in time for a crossing train when 20 < CUTOFF < 40. 

• Nuclear reactor controller. There are m rod-processes and one controller process. Each process has a 
clock with rate in [1, 1]. A rod just-moved out of the heavy water must stay out of water for at least 
T (a static parameter) time units. The timing constants used in the benchmarks are 58/10, 59/10, 16, 
and 161/10. The controller may miss the timing-constraints for the rods if — T < —(109m — 29)/5. 

• CSMA/CD. This is modified from [27]. The two timing constants A and B, set to 26 and 52 
respectively, are now treated as static parameters to be analyzed. We do require that B > 52. 
Basically, this is the ethernet bus arbitration protocol with the idea of collision-and-retry. The 
biggest timing constant used is 808. We want to verify that mutual exclusion after bus-contending 
period can be violated if A > A B > 52 A B < 808 A B < 2A. 

In our experiment, we compare performance in both forward and backward reachability analyses. The 
performance data of HyTech 2.4.5 and red 5.0 with dictionary ordering (no PSPSC), coefficient ordering 
(no PSPSC), magnitude ordering (PSPSC), and coefficient ordering with PSPSC is reported in table 3 (for 
backward analysis) and table 4 (for forward analysis). We stop experimenting with higher concurrency 
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benchmarks 


L.U11L. Lll 1 l^llL. V 


2 4 5 
(backward) 


red 5.0 (backward) 


dictionary | cocfFicicnt | magnitude 


coefficient 


no PSPSC 


DCDCP 


Fischer's 
mutual 
exclusion 
(m 


2 procs 


0.23s 


0.10s/17k 


0.11s/17k 


0.11s/17k 


0.07s/16k 


3 procs 


2.40s 


"1 o foil 

1.83s/81k 


1.75s/ 74k 


1.23s/59k 


O 7Tlc. I A Ah 
U. f US/ 44K 


A T""i v / \ r* o 


28.04s 


zu. zys / ozuk 


ZO.oDS / ZOyK 


IZ.OoS / Z1DK 


D. 14S / LOOK 


5 procs 


D /A/1 

W / Ivi 


278.8s/1420k 


354.1s/1149k 


162.0s/1034k 


Q1 'ifio / A 7 Ah 
Ol.OOS / 4 / 4K 


6 procs 


D / A/1 

W / Ivi 


2846s/5848k 


9923s /8796k 


1485s/4000k 


1 fi£ fit; /1 1 7Dk 


general 
railroad 
crossing 


2 trains 


O/M 


0.79s/103k 


0.68s/101k 


0.68s/101k 


n /n 1 l 
U. ( DS/y4K 


3 trains 


O / A/1 


11.48s/806k 


8.85s/616k 


8.84s/616k 


11/1 Q<n / e; qoi, 
11.4oS/ooUK 


4 trains 


D / A/1 

w / Ivi 


Z4o.DS/OU46K 


184. ys/ 4z4yk 


loD. IS/ 4Z4yK 


zoz.os/ ZoZUk 


5 trains 


O /M 


6095s/37093k 


4883s/25841k 


4900s/25841k 


ODZ f S/ lyZo4K 


reactor 
(m rods) 


2 rods 


0.056s 


0.08s/19k 


0.07s/19k 


0.06s/19k 


0.05s/15k 


3 rods 


0.33s 


0.41s/51k 


0.38s/52k 


0.37s/52k 


O OOc- / /I 1 \r 

U.ZZS f 41K 


4 rods 


2.61s 


3.10s/187k 


z.69s/186k 


z.71s/186k 


1.4ZS / 100K 


— 1 r °!j S — 


31 29s 


41.47s/104zk 


37.03s/ 1039k 


36.89s/ 1039k 


1 Q fi7<n /QC /ll, 
15. /S/oo4K 


o rods 


D^t I .OS 


951.5s/8228k 


866.9s/8191k 


839.3s/8191k 


4ul.oS/uy4lK 


CSMA/CD 
(m senders) 


2 senders 


O/M 


0.98s/42k* 


1.47s/125k 


0.57s/34k 


0.56s/33k 


3 senders 


O/M 


O/M 


5076s/ 2407k* 


121.5s/807k 


0.66s/105k 


4 senders 


O/M 


O/M 


O/M 


O/M 


2.47s/378k 


5 senders 


O/M 


O/M 


O/M 


O/M 


9.77s/1192k 


6 senders 


O/M 


O/M 


O/M 


O/M 


40.58s/3513k 



data collected on a Pentium 4M 1.6GHz with 256MB memory running LINUX; 
s: seconds; k: kilobytes of memory in data-structure; O/M: Out of memory; 

Table 3: Comparison in backward analysis with HyTech w.r.t. number of processes 



benchmarks 


concurrency 


HyTech 
2.4.5 

(forward' 


red 5.0 (forward) 


dictionary | coefficient | magnitude 


coefficient 


no PSPSC 


PSPSC 


Fischer's 
mutual 


2 procs 


0.34s 


0.10s/20k 


0.10s/20k 


0.10s/19k 


0.08s/18k 


3 procs 


37.89s 


O/M 


22.10s/561k 


19.18s/654k 


5.59s/538k 


general 
railroad 


2 trains 


3.29s 


2.29s/192k 


1.41s/95k 


1.43s/95k 


0.44s/84k 


3 trains 


O/M 


O/M 


O/M 


O/M 


6.35s/418k 


reactor 


2 rods 


O/M 


O/M 


O/M 


O/M 


O/M 


CSMA/CD 
(m senders) 


2 senders 


0.19s 


0.19s/29k 


0.17s/29k 


0.17s/29k 


0.25s/33k 


3 senders 


2.63s 


1.81s/102k 


1.64s/101k 


1.62s/101k 


2.61s/106k 


4 senders 


68.75s 


20.07s/370k 


17.49s/378k 


17.52s/378k 


27.03s/378k 


5 senders 


O/M 


268.0s/1905k 


240.3s/1906k 


242.2s/1906k 


331.9s/1910k 


6 senders 


O/M 


3889s/11725k 


3123s/11525k 


3155s/11525k 


4163s/11552k 



data collected on a Pentium 4M 1.6GHz with 256MB memory running LINUX; 
s: seconds; k: kilobytes of memory in data-structure; O/M: Out of memory; 

Table 4: Comparison in forward analysis with HyTech w.r.t. number of processes 

when we feel that too much time (like more than 1 hour) or too much memory (20MB) has been consumed 
in early fixpoint iterations. The experiment, although not extensive, does show signs that HRD-technology 
(with or without PSPSC) can compete with the technology used in HyTech 2.4.5. For all the benchmarks, 
HRD-technology demonstrates better scalability w.r.t. concurrency complexity. We believe that the data- 
sharing capability of HRD, when properly programmed, is the main reason for the performance advantage 
in the experiment. 

Finally, PSPSC cuts down the time and memory needed for parametric analysis. Especially, in forward 
analysis of the general railroad benchmark with three trains, without PSPSC, the state-space exploration 
fails to converge. This shows very good promise of this technique. 

12.2 Comparison with TReX 1.3 

Another famous tool for the verification of hybrid systems is TReX [1,3], which supports the verification 
of systems with clocks, static parameters, and lossy channels. As mentioned in section 2, the time-progress 
weakest precondition (or strongest postcondition in forward analysis) calculation algorithm in red 5.0 is 
more complex than the one in TReX. And TReX now mainly runs in forward analysis. And TReX also 
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benchmarks 


V_. Wilt., HI I V 


Forward 


Backward 


1 JxCA 

1 Q 
l.O 


red 5.0 


1 rveA 

1 Q 
l.O 


red 5.0 


magnitude 


LOt!ii . n^-T or uL 


magnitude 


LOtill . n^ir Or uL 


Fischer's 
mutual 
exclusion 
(m 


2 procs 


1.12s 


n 07.- /1 71.- 

U.U f S/ 1 / K 


U.U / S/ lOK 


8.96s 


n nor. /i qi, 
U.Uos/ loK 


U.U4S/ lok 


3 procs 


r\ /-\a 


l.OOS/ lo 1 K 


n 7s2 t -. / 701. 

U. foS/ ( yK 


\ A 1 Q7c 

>4iy < s 


U.DDS / 4oK 


U.4ys/4ok 


A nrnrc 
'-t ULn 


O / A/1 
VJ/ 1VJ 


1 Q7 Oc- /Q71 /IL 

iy / .ys / z i i4tK 


1 OOc- / t^QOL- 

iD.yzs/ooyK 


M /A 


Qic /1 sol 

D.olS/ loUK 


Q c;Cc /i c;ci- 
O.OOS/ lOOK 


5 procs 


O /A/1 


> loUUS 


(OZ. IS/ ozo4k 


1ST / A 


oy.z / s/y4oK 


Z4. / IS/DOOK 


6 procs 


/T\/i 
U/iVl 


NT / A 


\1 QOflc 

>loUU£ 


1V7 / A 
IN / A 


□0 / .1S/4o41K 


1 7fl /Q7GSL 
1 / U.oS/ Z ( yoK 


reactor 
(m rods) 


2 rods 


O/M 


O/M 


O/M 


N/A 


0.06s/19k 


0.05s/15k 


3 rods 


O/M 


O/M 


O/M 


N/A 


0.37s/52k 


0.22s/41k 


4 rods 


O/M 


O/M 


O/M 


N/A 


2.71s/186k 


1.42s/155k 


5 rods 


O/M 


O/M 


O/M 


N/A 


36.89s/1039k 


18.67s/884k 


6 rods 


O/M 


O/M 


O/M 


N/A 


839.3s/8191k 


461.8s/6941k 



Data for TReX (backward analysis) is collected on a Pentium III lGHz/900MB running Linux with CPU 

time normalized with factor 1/1.6. 
Data for red and for TReX (forward analysis) is collected on a Pentium 4M 1.6GHz/256MB running 

LINUX. 

s: seconds; k: kilobytes of memory in data-structure; 
O/M: Out of memory; N/A: not available; 



Table 5: Performance comparison with TReX w.r.t. number of processes 



may have tuned its performance for systems with lossy channels. Thus it can be difficult to compare 
the performance of TReX with red 5.0 directly. Anyway, we still tried hard and used one week to learn 
the input language of TReX and to analyze two benchmarks. The first is Fischer's protocol with all 
clocks in the uniform rate of 1. The second is the Nuclear Reactor Controller. The performance data is 
shown in table 5 for both forward and backward analysis. Two additional options of red 5.0 were chosen: 
coefficient evaluation ordering with PSPSC and magnitude evaluation ordering without. At this moment, 
since we do not have the reduce library, which is not free, in TReX for backward analysis, TReX team has 
kindly collected TReX's performance in backward analysis for us. Although the data set is still small and 
incomplete, but we feel that the HRD-technology shows a lot of promise in the table. We believe this can 
largely be attributed to the data-sharing capability of BDD-like data-structures. 



13 Summary 

This work is a first step toward using BDD-technology for the verification of LHAs. Although the initial 
experiment data shows good promise, we feel that there are still many issues worthy of further research 
to check the pros and cons of HRD-technology. Especially, we have to admit that we have not developed 
algorithms to eliminate general redundant constraints in HRDs. Our present implementation eliminates 
redundant LH-constraints that can be deduced by four peer LH-constraints along the same paths. We 
also require that the LH-expression of the redundant LH-constraint must not precede the LH-expressions 
of these four peer LH-constraints. Although our current implementation does perform well against the 
benchmarks, we still hope that there is a better way to check redundancy. 

Also, subsumption is another challenge. Straightforward implementation may use the complement 
of the current reachable state-space to filter those newly constructed weakest preconditions. Since the 
HRD of the current reachable state-space can be huge, its complement is very expensive to construct and 
maintain. 
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